GeRgraphy

16 March, 2016

GIS

Geographic Information Systems

GIS

Geographic Information Systems

incorporating
storing
manipulating
analyzing
displaying…

Spatial Data

What is spatial data?

Nonspatial data has no location information

nonspatial = data.frame(
  id=c(1,2,3,4),
  data=rnorm(4)
)
print(nonspatial)

##   id       data
## 1  1 -0.1140578
## 2  2  1.0914729
## 3  3 -0.4780712
## 4  4  0.1558649

What is spatial data?

Spatial data has location information

The simplest spatial data are points on a map

spatial = data.frame(
  id=c(1,2,3,4),
  data=rnorm(4),
  x=runif(4,-180,180),
  y=runif(4,-90,90)
)
print(spatial)

##   id        data           x         y
## 1  1 -1.06768167    8.198913  45.85735
## 2  2 -0.36974434  103.754462 -84.13661
## 3  3  0.03500162 -124.050690  39.74547
## 4  4  0.41151930   43.201727 -40.59028

What is spatial data?

Which we can convert to explicitly spatial data using the sp package. Most GIS packages in R store data as sp classes.

library(sp)

What is spatial data?

The sp package has a method called coordinates that converts points to an sp class.

coordinates(spatial) = ~ x + y
class(spatial)

## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"

plot(spatial, axes=T)

What is spatial data?

Spatial data also needs a reference system or "projection" that allows us to represent spatial features on a map. Projections can be thought of as simply a coordinate system with an origin that is relative to a known point in space.

This is a whole field of mathematically intensive study termed "geodesy"

Much of the field of geodesy is jam-packed in the rgdal package, which is a wrapper for the Geospatial Data Abstration Library

library(rgdal)

What is spatial data?

rgdal includes a comprehensive list of projections that are typically represented as a string of parameters.

The most common is our standard latitude/longitude system, where the coordinates are angular and the origin is the equator directly south of Greenwich, England. The simplest projection string to denote this projections is:

"+proj=longlat"

To define the projection for spatial, we write to its proj4string slot:

proj4string(spatial) = "+proj=longlat"

Projections are a necessary evil for GIS users (to be continued)

What is spatial data?

With a projection associated with our spatial data, we can now relate it to other spatial data. In other words, let's make a map!

library(leaflet)

m = leaflet(data=spatial) %>%
  addTiles() %>%
  addMarkers()
m

What is spatial data?

Arguably, the most common use of spatial data, and therefore GIS, are in data exploration, simply taking a look at where the phenomenon occur. This is useful in the sense that we can see obvious trends in the data, giving clues for further analysis, or merely checking our work, seeing if we have made mistakes. Additionally, spatial data is used for different kinds of analysis, leveraging the location data to make inferences. This analysis could take the form of selecting all wells that fall within certain soil types or in the construction of error terms in construction a statistical model.

Maybe delete this paragraph… The most common software used in GIS are ArcGIS, produced by ESRI and QGIS, an open source analog, but there are many other software packages for GIS, both propietary and open-source. Any familiarity with the most popular GISs will show you that they are GUI-based and therefore much different from the look and feel of R, and, while R is not necessarily a replacement for the ArcGIS or QGIS, the analysis and plotting capability of R makes it a powerful GIS, and when it comes to model construction, R, and no suprise here, blows the common GIS platforms out of the water.

Spatial data types

Two main types: vector and raster

Vector = Polygons
Raster = Grid

Vector = Discrete
Raster = Continuous

Vector = Illustrator/Inkscape
Raster = Photoshop/GIMP

Maybe break this apart and put this into the respective intros of the datatypes below? Spatial data are divided up into two kinds of data, vector data and raster data. Vector data consist of the location data being defined by points. There are three main kinds of vector data: points, lines and polygons. Point data, as might be guessed, are data where each observation is defined by one coordinate pair, X and Y. In line data each observation is defined by a series of coordinate pairs that define a one dimensional path in space. Polygons are data where each observation is defined by coordinate pairs that trace out a two-dimensional area. Common uses of vector data are for phenomenon that have discrete boundaries and locations, e.g., wells and bird sitings for points, rivers and streets (at certain scales) for lines, and county borders and soil types for polygon data. Raster data are constructed as grids that are laid out over an area. Each cell of the grid contains some information about the phenomenon, but each cell will not be explicitly given a coordinate pair, rather every cell has the same length and width and so only the certain corners of the data need to be given coordinate pairs. Precipitation and temperature data are often portrayed as raster data because of the way they vary continuously over the space. Examples of the creation, input, and output of these data types are given below.

Note to power users: The usual definitions of these data describe vector data as being for data with discrete boundaries and raster data for that which continuously varies. While this is true in general, it is essentially a matter of convenience that determines which type is best. For example, you could represent a dataset of bird sitings in a raster format, but because the grid is so sparse (NA or 'no siting' for nearly all the cells) it saves memory (and makes more inuitive sense) to store the data as points. Conversely, rainfall data can be stored as thousands of points, each marking the estimated rainfall at that location, but storing the data in this format is rather cumbersome if the area is large.

This presentation is meant to give an introduction to handling spatial data in R. It's intended audience are those with at least a basic knowledge of R and an interest in, though not necessarily knowledge of, working with spatial data. We will cover basic handling and manipulation of vector data and raster data, and then give a flyover of how coordinate reference systems work and how to work with them in R.

Vector Data

Intro

Vector I/O

library(rgdal)

Vector I/O

soils = readOGR(dsn="data", layer="soilsData")

## OGR data source with driver: ESRI Shapefile 
## Source: "data", layer: "soilsData"
## with 75 features
## It has 27 fields

writeOGR(
  soils,
  "data",
  "soilsData_out",
  driver="ESRI Shapefile"
)

Vector data structure

class(soils)

## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"

slotNames(soils)

## [1] "data"        "polygons"    "plotOrder"   "bbox"        "proj4string"

Vector data structure

Take a peak here at the top of our attribute table

str(soils@data[,1:5])

## 'data.frame':    75 obs. of  5 variables:
##  $ mukey  : Factor w/ 25 levels "2774742","2774772",..: 12 19 12 6 9 5 7 24 25 13 ...
##  $ muarcrs: Factor w/ 75 levels "0.40538405","0.90105194",..: 25 62 70 19 3 29 20 7 5 27 ...
##  $ Sand1  : num  21.5 11 21.5 12 29.5 ...
##  $ Sand2  : num  35.29 8.34 35.29 9.02 38.22 ...
##  $ Sand3  : num  45.09 7.56 45.09 16.59 64.52 ...

*Make sure not to call head() or str() on a spatial object, as this spits out the first six features and all their attributes

Vector data structure

str(soils@data[,1:10])

## 'data.frame':    75 obs. of  10 variables:
##  $ mukey  : Factor w/ 25 levels "2774742","2774772",..: 12 19 12 6 9 5 7 24 25 13 ...
##  $ muarcrs: Factor w/ 75 levels "0.40538405","0.90105194",..: 25 62 70 19 3 29 20 7 5 27 ...
##  $ Sand1  : num  21.5 11 21.5 12 29.5 ...
##  $ Sand2  : num  35.29 8.34 35.29 9.02 38.22 ...
##  $ Sand3  : num  45.09 7.56 45.09 16.59 64.52 ...
##  $ Sand4  : num  48.6 30.6 48.6 36.7 33.2 ...
##  $ Sand5  : num  0 31.1 0 27.1 57.2 ...
##  $ Silt1  : num  45.8 65.2 45.8 68.7 54.5 ...
##  $ Silt2  : num  37.3 64.9 37.3 60.3 43.5 ...
##  $ Silt3  : num  36.7 64.1 36.7 34 21.1 ...

Vector data structure

str(soils@polygons[1])

## List of 1
##  $ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. ..@ Polygons :List of 1
##   .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
##   .. .. .. .. ..@ labpt  : num [1:2] -90.1 43.1
##   .. .. .. .. ..@ area   : num 1.19e-06
##   .. .. .. .. ..@ hole   : logi FALSE
##   .. .. .. .. ..@ ringDir: int 1
##   .. .. .. .. ..@ coords : num [1:21, 1:2] -90.1 -90.1 -90.1 -90.1 -90.1 ...
##   .. ..@ plotOrder: int 1
##   .. ..@ labpt    : num [1:2] -90.1 43.1
##   .. ..@ ID       : chr "0"
##   .. ..@ area     : num 1.19e-06

Vector data structure

Number of features

length(soils)

## [1] 75

Vector data structure

Index the first feature with a slice

poly_1 = soils[1,]

Plotting

plot(
  soils,
  main="Soils from Western WI",
  col=rainbow(5)
)

Working with Vector Data

Raster Data

Intro

A raster grid is rectangular.

Grid is another word for matrix.

Grid is another word for image.

A GIS raster grid is a matrix/image with an associated location and projection.

Intro

At a minimum, a GIS raster grid contains:

matrix of values
projection
reference point, often (x,y) of the lower-left corner
cellsize

Raster I/O

The rgdal rgdal packages is primarily for I/O and projecting GIS data

library(rgdal)

The raster package does everything rgdal does, but it includes lots of additional functionality.

library(raster)

Raster I/O

elev = readGDAL("data/dem_wi.tif")
writeGDAL(elev, "data/dem_wi_out.tif")

elev = raster("data/dem_wi.tif")
writeRaster(elev, "data/dem_wi_out.tif")

Raster data structure

The raster object elev has all the necessary pieces of spatial information:

elev

## class       : RasterLayer 
## dimensions  : 284, 387, 109908  (nrow, ncol, ncell)
## resolution  : 0.01666667, 0.01666667  (x, y)
## extent      : -93.03262, -86.58262, 42.3949, 47.12823  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## data source : C:\Users\ruesca\Documents\GeRgraphyPresentation\data\dem_wi.tif 
## names       : dem_wi 
## values      : 175, 565.4104  (min, max)

Raster data structure

Which means we can make a map!

m = leaflet() %>%
  addTiles() %>%
  addRasterImage(elev, opacity=0.5)
m

Raster data structure

Raster analysis

Remember that rasters are just matrices!

Therefore, most matrix operations can be applied to rasters. For example:

plot(
  elev > 400,
  col=c("red", "blue")
)

Raster analysis

Rasters can be easily converted to matrices to do more complex work.

lat_grad = apply(
  as.matrix(elev),
  1,
  mean,
  na.rm=T
)
plot(lat_grad, type="l")

Raster overlay

Most raster analysis ultimately executes some sort of overlay.

The issue:

To overlay two or more rasters, their projections, extents, and cellsizes must align perfectly.

This can be a difficult task.

Raster overlay

coordinate systems

What is the highest point in each county?

# Pseudo-code
1. Read in elevation data (raster grid)
2. Read in county boundary data (polygons)
3. Convert counties to raster grid that aligns with elevation grid
4. Find maximum elevation gridcell within each county

Raster overlay

coordinate systems

counties = readOGR("data", "WI_Counties")

## OGR data source with driver: ESRI Shapefile 
## Source: "data", layer: "WI_Counties"
## with 72 features
## It has 7 fields

elev

## class       : RasterLayer 
## dimensions  : 284, 387, 109908  (nrow, ncol, ncell)
## resolution  : 0.01666667, 0.01666667  (x, y)
## extent      : -93.03262, -86.58262, 42.3949, 47.12823  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## data source : C:\Users\ruesca\Documents\GeRgraphyPresentation\data\dem_wi.tif 
## names       : dem_wi 
## values      : 175, 565.4104  (min, max)

Raster overlay

coordinate systems

proj4string(counties)

## [1] "+proj=tmerc +lat_0=0 +lon_0=-90 +k=0.9996 +x_0=520000 +y_0=-4480000 +ellps=GRS80 +units=m +no_defs"

proj4string(elev)

## [1] "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"

Raster overlay

coordinate systems

extent(counties)

## class       : Extent 
## xmin        : 294839 
## xmax        : 770036.4 
## ymin        : 225108.8 
## ymax        : 734398.4

extent(elev)

## class       : Extent 
## xmin        : -93.03262 
## xmax        : -86.58262 
## ymin        : 42.3949 
## ymax        : 47.12823

Raster overlay

coordinate systems

cty_grid = rasterize(counties, elev, field="COUNTY_FIP")
summary(cty_grid)

##          layer
## Min.        NA
## 1st Qu.     NA
## Median      NA
## 3rd Qu.     NA
## Max.        NA
## NA's    109908

Raster overlay

coordinate systems

prj = proj4string(elev)
cty_prj = spTransform(counties, prj)

To do this, we use the spTransform function in the sp package.

Raster overlay

coordinate systems

extent(cty_prj)

## class       : Extent 
## xmin        : -92.88924 
## xmax        : -86.8048 
## ymin        : 42.49197 
## ymax        : 47.08077

extent(elev)

## class       : Extent 
## xmin        : -93.03262 
## xmax        : -86.58262 
## ymin        : 42.3949 
## ymax        : 47.12823

Raster overlay

coordinate systems

plot(elev)
plot(cty_prj, add=TRUE)

Raster overlay

coordinate systems

cty_grid = rasterize(cty_prj, elev, field="COUNTY_FIP")
summary(cty_grid)

##         layer
## Min.        1
## 1st Qu.    18
## Median     36
## 3rd Qu.    55
## Max.       72
## NA's    50557

Raster overlay

coordinate systems

extent(cty_grid)

## class       : Extent 
## xmin        : -93.03262 
## xmax        : -86.58262 
## ymin        : 42.3949 
## ymax        : 47.12823

extent(elev)

## class       : Extent 
## xmin        : -93.03262 
## xmax        : -86.58262 
## ymin        : 42.3949 
## ymax        : 47.12823

Raster overlay

coordinate systems

library(dplyr)
ovly = data.frame(
  elev=getValues(elev),
  cty=getValues(cty_grid)
)

hi_pt = ovly %>%
  group_by(cty) %>%
  mutate(
    elev = (elev == max(elev, na.rm=T)) * elev
  ) %>%
  ungroup()

elev = setValues(elev, hi_pt[["elev"]])
elev[elev == 0] = NA
hi_pt_sp = rasterToPoints(elev, spatial=T)